Anomalous shell effect in the transition from a circular to a triangular billiard 
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We apply periodic orbit theory to a two-dimensional non-integrable billiard system whose bound- 
ary is varied smoothly from a circular to an equilateral triangular shape. Although the classical 
dynamics becomes chaotic with increasing triangular deformation, it exhibits an astonishingly pro- 
nounced shell effect on its way through the shape transition. A semiclassical analysis reveals that this 
shell effect emerges from a codimension-two bifurcation of the triangular periodic orbit. Gutzwiller's 
semiclassical trace formula, using a global uniform approximation for the bifurcation of the trian- 
gular orbit and including the contributions of the other isolated orbits, describes very well the 
coarse-grained quantum-mechanical level density of this system. We also discuss the role of discrete 
symmetry for the large shell effect obtained here. 

PACS numbers: 05.45.Mt, 03.65.Sq 



I. INTRODUCTION 

The periodic orbit theory (POT) is a very useful tool to 
study shell structure in single-particle energy spectra. In 
POT, the quantum-mechanical level density is semiclassi- 
cally approximated in terms of the periodic orbits of the 
corresponding classical system. For systems with only 
isolated orbits, Gutzwiller derived the so-called "trace 
formula" |l[ which is particularly successful for chaotic 
systems. The POT for 3D cavities was developed in |2j. 
In integrable systems, semiclassical trace formulae can 
be derived from torus quantization 0,0). However, most 
physical systems lie between the above two extreme sit- 
uations, i.e., they exhibit mixed phase-space dynamics 
in which both regular and chaotic motion coexists on the 
same energy shell. For systems with various types of con- 
tinuous symmetries, trace formulae have been derived in 
0) & @ • For an introductory text book on semiclassical 
physics and applications of the POT to various physical 
systems, we refer to Q- 

In both integrable and mixed systems, bifurcations of 
periodic orbits can significantly influence the shell struc- 
ture. The above-mentioned semiclassical trace formulae, 
which are all based on the stationary-phase approxima- 
tion (SPA), diverge when bifurcations of periodic orbits 
occur under variations of energy or potential parame- 
ters (e.g., deformations). In the SPA, the classical ac- 
tion integral is expanded around its stationary points 
(corresponding to periodic orbits) up to quadratic terms 
and the trace integral is evaluated in terms of Gauss- 
Fresnel integrals. At bifurcation points, the determinant 
of the coefficient matrix of the quadratic terms becomes 
zero and the Gauss-Fresnel integral becomes singular. In 
order to obtain a finite semiclassical level density near 
bifurcation points, higher-order expansion terms of the 
action integral are needed, which are most conveniently 
found from the normal forms appropriate to the types of 
bifurcation under consideration [g, |9|, [10| . 

Quantum billiards (and their 3D versions, cavities), be- 
sides being quite useful toy models to study POT, reflect 



important features of finite physical quantum systems 
such as quantum dots, metallic clusters, and atomic nu- 
clei. E.g., non-integrable 3D cavities with realistic shapes 
appropriate for fission barriers of actinide nuclei have 
been used in POT to explain the onset of the mass asym- 
metry of the nascent fission fragments [ll[ . On the other 
hand, many integrable billiard systems are well known 
and fully understood semiclassically, e.g., circular, equi- 
lateral triangular, square, and elliptic billiards (cf. (7j), 
and may be used as simple models for physical systems. 
In these, bifurcations of short periodic orbits may lead 
to a considerable enhancement of shell effects. E.g., in 
the elliptic billiard, the short diametric orbit undergoes 
successive bifurcations with increasin g de formation, and 
new periodic-orbit families emerge 

[iHlGl- The same 
type of bifurcations occur for equatorial orbits in the 
3D spheroidal cavity [3] and provide a schematic ex- 
planation of nuclear sup er-deformed and hyper-deformed 
shell structure [HI, llfjj ]. In these studies, it was shown 
that a system may turn strongly chaotic by adding small 
reflection-asymmetric (e.g., octupole) deformations (l7| . 



In this paper, we apply the POT to a two-dimensional 
non-integrable billiard system whose boundary is contin- 
uously varied from a circular to an equilateral triangular 
shape. This study is initiated to explore a possible quan- 
tum dot system. Many studies have been undertaken for 
quantum dots, but this type of deformation has not been 
studied before. Another important aim is to investigate 
the role of discrete symmetries. As will be discussed, 
the present model possesses the discrete C^ v point-group 
symmetry. Recently, several mean-field studies of nuclei 
have suggested the possible existence of low-lying states 
with exotic shapes with point-group symmetries, such as 
tetrahedral and octahedral deformations [lSj. In order 
for such shapes to be stabilized, rather large quantum 
shell effects in the single-particle spectra are necessary. 
We will discuss the role of discrete symmetries in estab- 
lishing strong shell effects. 



2 



fyi=0+ 0- 1+ 1- 




FIG. 1: Shapes of the boundary given by Eq. (|2.2p for several 
values of deformation parameter a. 



II. SHELL STRUCTURE AND LEVEL 
STATISTICS 

A. The model system 

We consider the two-dimensional billiard system 

H = f + V(r), V(r) = {° r < Jfjj (2.1) 
2m {°° r > R{9) v ' 

in polar coordinates r = (r, 0) , whereby the boundary 
shape R(9) is parameterized implicitly by the equation 

R 2 + ^^^_ cos{3d)=R 2 j ee[Q,2n). (2.2) 

Figure [1] shows the shape of the boundary for several val- 
ues of a. This system possesses the discrete symmetries 
of the point group Ca v , consisting of ±27r/3 rotations 
about the origin and reflections with respect to the three 
axes through the origin with 9 = 0, ±7r/3. The deforma- 
tion parameter a = yields a circular shape and a = 1 
an equilateral triangle. The system is integrable in these 
two limits, but non-integrable in between. 

For the calculation of the quantum spectrum of this 
system, it is useful to decompose the wave functions into 
free circular waves: 

oo 

Mr,0)= J2 c m J\ m \{kr)e lm \ (2.3) 

rn — — oo 

where k = \/2me/h is the wave number and e the energy. 
Taking the C^ v symmetry into account, we can classify 
the eigenstates according to the eigenvalues of the sym- 
metry operators 

TfyO*) = e 27riA/3^(AM) ) ( 24a ) 

■p^) = (_ 1 )P^(Am) j (2.4b) 



08.143 07.733 



09.338 09.181 07.979 

FIG. 2: Wave functions of some lowest eigenstates for a = 0.4. 
0+, 0- and 1± are states of type (f23a)l . (|2.5b|) and (f23cT ). 
respectively. Their k values are also indicated. 

where 1Z and V represent rotation by 27t/3 around the 
origin and reflection with respect to the x axis, respec- 
tively. We have simultaneous eigenstates of 1Z and V for 
A = states, and then we can classify the eigenstates 
into four sets, which correspond to the irreducible repre- 
sentations ( "irreps" ) of the C 3v point group [HI : 

oo 

V4° +) (r, 0)=Y1 c ™ +) J 3m(fcr) cos(3m0) , (2.5a) 

rn— 

oo 

^ 0-) (r> 9) = J2 cgP-Wfcr) sm(3m9) , (2.5b) 

771—1 

OC 

^ ±1} M)= E c^J l3m±ll (kry^ e . (2.5c) 

m— — oo 

V>1 are the complex conjugates of ipu~ , and the two 
form degenerate pairs of states. (The point group C 3v 
has two 1-dimcnsional irreps and one 2-dimensional irrep. 
The states (0±) correspond to the 1-dimensional irreps, 
and the degenerate pairs of states (±1) correspond to 
the 2-dimensional irrep.) Taking linear combinations of 
these, one finds the following alternative real expressions 
for the states (I2.5c|) : 

m— — oo ** J 

EM) 

which are not eigenstates of the operator 1Z, but of the 
operator V . The eigenvalue spectrum {k n } and the co- 
efficients c m are determined by the Dirichlet boundary 
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described by the Hamiltonian (|2.ip . We split the energy 
into a smooth and an oscillating part 



E(N) 



N/2 
n=l 



E(N) + 5E(N) , (N even) (2.6) 



whereby the factor 2 accounts for the spin s = 1/2. The 
smooth part E(N) is equivalently given by a Strutin- 
sky averaging [2(|, the extended Thomas- Fermi model, 
or the Weyl expansion (cf. 0], Ch. 4), while the shell- 
correction energy SE(N) reflects the quantum effects; it 
is dominated by the shortest periodic orbits of the clas- 
sical system as demonstrated below for the case of the 
level density. Large gaps in the spectrum lead to large 
amplitudes of 5E(N). This is demonstrated in Figure 
[4] where we present SE(N), scaled by a factor y/~N, for 
three values of a. We note that the oscillatory pattern 
for a = 0.5 is quite regular and on the average has a 
much larger amplitude than in the integrable limits. 



0.4 0.6 
a 



B. Level statistics 



FIG. 3: Lowest part of energy spectrum (in units of ft 2 /mR 2 .), 
plotted against the deformation parameter a. 




FIG. 4: Scaled shell-correction energy 8E(N) (in units of 
h 2 /mRo) versus y/N. 



condition ij}k n (R,0) — 0. We show some eigenfunctions 
in Fig. [2j The states belonging to the sets (|2.5aP and 
(|2.5bp have symmetry under the rotation 1Z. The former 
are even under the reflection V, while the latter are odd. 

Figure [3] shows the lower part of the energy spectrum 
{ e n} — {?i 2 ^/2m}, plotted against the deformation pa- 
rameter a. We note that the levels bunch into small en- 
ergy intervals in the region a ~ 0.5 — 0.7, forming large 
gaps in the spectrum, ft is, in fact, quite surprising to 
realize that the shell gaps in this non-integrable region 
are much larger than in the two integrable limits a — 
and a = 1. The gaps cause large shell effects in the 
total energy of a system of N noninteracting fermions 



Nearest-neighbor spacing (NNS) distributions are com- 
monly used to identify signatures of chaos in a quantum 
system. Generically, classically chaotic systems exhibit 
level repulsion and the NNS distributions are of Wigner 
type, while regular systems typically have degeneracies 
and the NNS distributions are Poisson like 2l|. To ex- 
tract these universal fluctuation properties, one has to 
use unfolded spectra whose mean level density is normal- 
ized to unity. Thus, for systems with discrete symme- 
tries, one has to study the NNS of the subsets of levels 
belonging to the irreps of the corresponding point group. 

The average total level density of a two-dimensional 
billiard is given by Weyl's asymptotic formula [22] 



9(e) 



2irk- 



■A + Oie- 1 ' 2 ), 



(2.7) 



in which the leading-order term is proportional to the 
area A surrounded by the boundary and does not depend 
on the energy. This means that for large energies e the 
mean level spacing A becomes asymptotically constant: 



A„ = i 
9 



mA 



(2- 



As discussed in the previous subsection, the quantum 
levels of our system fall into four sets {e^} with n = 0± 
and ±1, corresponding to the eigenstates (|2.5b .b) and 
(|2.5c[) or (|2.5c[ L respectively. The numbers of levels in 
these sets have the relative ratios 

iV(°+) :JV<°-> : JV<+ 1 ) rJV^ =1:1:2:2, 

and hence the mean level spacing in each set is given by 

A( 0± )=6A , A( ±1 )=3A . (2.9) 
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a=0.00 



a=0.50 




FIG. 5: Nearest-neighbor spacing distribution P(s) for four 
values of the deformation parameter a. The lowest 600 levels 
(i.e., 100 and 200 levels of the subsets 0± and ±1, respec- 
tively) were used to obtain the statistics. Solid and dashed 
lines show Poisson and Wigner distributions, respectively. 
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FIG. 6: Brody parameter to for nearest-neighbor spacing dis- 
tribution, plotted as a function of a. 



The unfolded level spacings are then obtained by 



„(«) 



'n+l 



AW 



(« = 0±,±1) 



(2.10) 



Figure [5] shows the NNS distributions P(s), averaged 
over all four sets. At a = 0, the system is integrable and 
the distribution is Poisson like as expected. At a = 0.3, 
the distribution changes into Wigner form. A similar sit- 
uation is also found at a = 0.7. At a = 0.5, however, the 
distribution deviates considerably from the Wigner form 
and becomes closer to a Poisson distribution. For mixed 
systems, the NNS distributions P(s) are often fitted by 
a Brody distribution [23| 



B(s,u) =/3(w + l)s^exp(-/3s" +1 ), (3 ■ 

(2.H) 

which interpolates between the Poisson (uj = 0) and 
Wigner (u> = 1) distributions. The Brody parameter 
li can then be taken as a measure for the chaoticity of 
NNS distribution. In Figure [5] we show u as obtained 
by fitting the P{s) distributions of Fig. [5] to (|2~TT|) . We 
clearly recognize two peaks around a ~ 0.3 and ~ 0.7, 



exhibiting near-chaoticity, separated by a deep minimum 
around a ~ 0.5 where the system appears to approach 
regularity As we will discuss below, this near-regularity 
is related to an approximate restoration of local dynam- 
ical symmetry due to a periodic-orbit bifurcation. 



III. FOURIER ANALYSIS AND CLASSICAL 
PERIODIC ORBITS 

In the POT, the quantum level density g(e) is approx- 
imated in terms of classical periodic orbits by the semi- 
classical trace formula [H, 0, S H @| 



5(e) = X^( e ~ e ") 



Ke) + ^A c (e)cos 



(3.1) 



where the first term, like for the energy in (|2.6p . repre- 
sents the smooth part, while the second term contains 
the quantum shell effects. In the latter, the sum is taken 
over all periodic orbits £ of the classical system (or the 
orbit families £ in systems with continuous symmetries 
[E 0, IE] ) , S{ is the action integral around £, and is the 
Maslov index [24], Hf| . The amplitude depends on the 
stability of the orbit £ (and, for an orbit family, on the 
phase-space volume covered by the family). For isolated 
orbits, the amplitude A^ was given by Gutzwiller [l|: 



At 



7f(e) 



irh ^/|dct[l -M ? (e) 



(3.2) 



where T^(e) = dS^(e)/de is the period of the orbit and 
Mjj(e) its stability matrix defined below. If an orbit has 
a discrete degeneracy / (i.e., if there exist / replicas with 
identical actions, stabilities and Maslov indices but dif- 
ferent orientations) due to discrete symmetries, it has to 
be included / times into the sum in (|3.1[) . This holds 
also for time- reversed rotational orbits. 

Transforming variables from energy e to wave number 
k and using the relation S% — hkL^ for billiards, where 
I/£ is the length of the orbit £, the trace formula becomes 



h 2 k 

rn 



9(e) 



s 



(3.3) 



Let us now consider the Fourier transform of the level 
density g(k) with respect to k: 



F{L) 



g{k)e~ lkL e-^ 2 / 2 dk. (3.4) 



The Gaussian damping factor is included for the trun- 
cation of the high-energy part of the spectrum. If we 
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insert Eq. (|3.3[) and ignore the k dependence of the am- 
plitude factors (which are constant for isolated orbits 
in billiards and cavities), we obtain 

F SC (L) =F (L)+7r^ e -^/ 2 A^ A (i-^). (3.5) 

6a(x) is a normalized Gaussian of width A, which turns 
into Dirac's delta function in the limit A — > 0. Eq. (|3.5p 
indicates that the Fourier transform F(L) is a function 
with successive peaks at the lengths of the periodic or- 
bits L = Li£, with heights proportional to the amplitudes 
. We can therefore extract information about the clas- 
sical periodic orbits from the Fourier transform of the 
quantum-mechanical level density: 



10 r \F(L)\ (arb. unit) 



k„ = \/2me n /h . (3.6) 



Figure shows the modulus of this Fourier transform 
of our system versus the deformation parameter a. At 
a = 0, where the periodic orbits form degenerate families 
corresponding to rational tori, we find the peaks at the 
well-known orbit lengths of the circular billiard 



L vw — 2vRq sin(w7r/u) 



(3.7) 



Here w and v are the winding number around the ori- 
gin and the number of vertices (v > 2w) of each orbit, 
respectively. For instance, L = 4i?o for the diametric 
orbit (w = 1, v = 2), L = 3\/3-Ro for the triangular or- 
bit (w = l,v = 3), L = 4^/2R for the square orbit 
(w = 1, v = 4), and so on. With increasing a, we observe 
a dramatic enhancement of the peak height correspond- 
ing to the triangular orbit, L ~ 5.2Rq (and its second 
repetition, L ~ 10.4), starting around a ~ 0.4 and cul- 
minating around a ~ 0.55. As we shall see below, this 
can be traced back to a codimension-two bifurcation of 
the triangular orbit. 

In fact, for a > all rational tori of the circular bil- 
liard are broken into pairs of stable and unstable isolated 
orbits. With increasing a, bifurcations of the stable or- 
bits occur and new periodic orbits emerge, which makes 
the phase-space increasingly chaotic. In the following, 
we first demonstrate this for the two shortest orbits and 
then focus on the triangular orbit. 

The stability of a periodic orbit is described by the 
stability matrix Mj, which is defined by the linearized 
Poincarc map around the periodic orbit 



Me 



%(0),p(0)) ; 



(3.8) 



where (q(t),p(t)) are the coordinates and momenta per- 
pendicular to the periodic orbit £ as functions of time t, 
and TJr is the period of the orbit. In a two-dimensional au- 
tonomous Hamiltonian system is a symplectic (2x2)- 
matrix, and the stability of a orbit is easily identified by 
looking at its trace TrM^. For elliptic (stable) orbits, 
the eigenvalues of are of the form (e lv ,e~ m ) with 




FIG. 7: Modulus |F qm (L)| of Fourier transform of the 
quantum-mechanical level density versus deformation param- 
eter a. 
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FIG. 8: Trace of stability matrix Tr Mf for the diameter orbits 
(2A) and the second repetition of the stable diameter orbit 
(2A 2 ), plotted versus deformation parameter a. 



real v > 0, and thus -2 < TrMj < 2. For direct- and 
inverse-hyperbolic (unstable) orbits, the eigenvalues are 
(e u , e~ u ) and (— e", — e _u ), respectively, with real u > 0, 
and hence TrAf 5 > 2 and TrM 5 < -2. Bifurcations of 
isolated orbits occur whenever TrM^ = +2. 

Figures [H] and [5] show TrMg(a) for the two shortest 
pairs of periodic orbits in our system. In Fig. [51 the sta- 
ble branch (2A) of the diametric orbit is seen to undergo 
a period-doubling pitchfork bifurcation at a — 0.166, 
where a symmetric wedge-shaped orbit (4V) emerges. 
The latter undergoes a further pitchfork bifurcation at 
a = 0.227, where a pair of asymmetric wedge-shaped 
orbits (4U) emerges. 

The bifurcation scenario of the triangular orbits is 
more complicated. In the upper panel of Fig. [51 where 
TrMj is plotted against \fa, the stable triangular or- 
bit (3A) is seen to touch the critical line TrMj = +2 
near ^fa ~ 0.63 (a ~ 0.4) and to remain stable on ei- 
ther side. A pair of new stable (3D) and unstable (3C) 
triangular orbits emerge from the touching point. A mag- 
nification of the situation around a ~ 0.4, shown in the 
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FIG. 9: Same as Fig. [8] but for the triangular orbits, in the 
upper panel plotted versus y/a. The lower panel shows a 
magnification around the bifurcation of the stable orbit 3A 
occurring near a ~ 0.407. 



lower panel against a, reveals that the scenario consists 
of two connected near- lying bifurcations, also called a 
"codimension-two bifurcation". At a = 0.400, a tangent 
(saddle-node) bifurcation occurs, at which the new orbits 
3D and 3C are born. Shortly after this bifurcation, the 
(old) stable orbit 3A and the new unstable orbit 3C en- 
counter in a touch-and-go bifurcation at a = 0.407 (see 
Appendix [Al for its analytic value). Note that the new 
3C and 3D orbits do not possess C^ v symmetry, in con- 
trast to the old 3 A orbit, so that they occur in degenerate 
triplets obtained by successive rotations about 2tt/3, i.e., 
these orbits have a discrete degeneracy of / = 3. 

Figure [TU] shows excerpts of Poincare surfaces of sec- 
tion (6, ip) in the relevant regions. Here is the polar 
angle of a reflection point of an orbit at the boundary, 
while ip represents the reflection angle measured from 
the normal to the boundary at the reflection point. The 
three upper panels illustrate the tangent bifurcation of 
the 3C and 3D orbits. For a = 0.390 (upper left), one 
sees one major regular island corresponding to the stable 
equilateral triangular orbit 3A. It contains its fixed point 
at (6>, ip) = (0, 7r/6), surrounded by quasi-tori (small ape- 
riodic perturbations of the 3A orbit). At the bifurcation 
point a ~ 0.400 (upper center), three cusps are formed 
by one of the surrounding quasi-tori in the island, and for 



a = 0.401 (upper right), the three stable fixed points of 
the new 3D orbits, surrounded by small regular islands, 
are seen to have emerged from each of the three cusps 
in the major island. The three saddles separating these 
three islands from the central island contain the unstable 
fixed points of the 3C orbits. The stable fixed point of 
the 3A orbit still persists at the original position at the 
center of the major island. The number / = 3 of the new 
stable and unstable fixed points in the island is due to 
the three-fold discrete degeneracy of the orbits 3D and 
3C mentioned above. 

The lower three panels in Fig. ITDl illustrate the touch- 
and-go bifurcation of the orbits 3A and 3C. At the bi- 
furcation point a = 0.407 (lower left), the three unstable 
fixed points of the 3C orbits have contracted into a star- 
like intersection of three quasi-tori, located at the central 
fixed point (9, ip) = (0, 7r/6) of the 3A orbit. The three 
nearby stable fixed points in the major island belong to 
the orbits 3D. For a = 0.420 (lower center), small sta- 
ble islands have formed again around the fixed points of 
the 3 A orbit, and nearby we recognize the three unsta- 
ble fixed points of the 3C orbit. For a — 0.500 (lower 
right), the central island of the 3A orbits have grown, 
the three islands of the stable 3D orbits have shrunk, 
and the three unstable fixed points of the 3C orbits are 
about to be buried in the increasing chaotic structure. 

The bifurcation of the triangular orbit 3A occurs at the 
deformation where we observed the onset of the large 
Fourier peak in Fig. [7) It is due to the appearance of 
the pair of six- fold degenerate new orbits 3D and 3C. 
In particular, the stable one (3D), besides the doubly- 
degenerate 3 A orbit, adds to the local regularity of the 
phase space observed in Figs. 151 171 and [TU1 It is therefore 
apparent that this bifurcation is responsible for the re- 
markable shell structure in the quantum spectrum seen 
in Fig. [21 Note that the maximum of the Fourier peak 
in Fig. [7] occurs at a ~ 0.55, i.e., above the bifurcation 
where the new orbits are born. The fact that the sig- 
nature of a bifurcation is strongest after the bifurcation 
point (on the side where the new orbits exist) is a rather 
general trend for shell effects which can also be under- 
stood in terms of semiclassical uniform approximations 
@> 03 (see the following section and Fig. [12] below) . It 
has also been observed recently in the level statistics of 
a Hamiltonian system with mixed phase space [26| . 



IV. SEMICLASSICAL ANALYSIS 



A. Semiclassical level density near bifurcations 



As mentioned above, the Gutzwiller trace formula ( 
with the amplitudes (|3.2[) for isolated orbits diverges at 
bifurcation points. This is due to the break-down of sta- 
tionary phase approximation, and higher-order terms in 
the expansion of the action integral must be included in 
the derivation of the trace formula. For this purpose, 
it is necessary [8[ to express the trace integral in phase 
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FIG. 10: Poincare surfaces of section (9, sin ip) in the fundamental domain of the C3 group, plotted for several values of a (given 
at the upper left of each panel) in the bifurcation region of the triangular orbit 3A. 9 is the polar angle of a reflection point 
and ip the reflection angle measured from the normal to the boundary; the set (9, sin ip) is approximately area preserving. 



space. After integrating over the coordinate and momen- 
tum parallel to an isolated orbit £ in a two-dimensional 
system, its semiclassical contribution to the oscillating 
part of the level density becomes 0, 



%(e) 



2ir 2 h 2 



x exp 



Re 



dq< 



, 1 dS 

n r oe 



d 2 s 



dpdq' 



-S(q',p) - -q'p- 



-v 6 



1/2 



(4.1) 



Here q' is the final coordinate and p the initial momentum 
perpendicular to the orbit £, the two forming a canonical 
pair of variables (q',p) to describe the "natural" Poincare 
surface of section (PSS) of the orbit on which (q',p) — 
(0, 0) is its fixed point. For generic bifurcations, n r is the 
repetition number of the primitive orbit. The function 
S(q',p) denotes the Legendre transform of S(q',q): 



S{q',p) = S(q',q) +qp, 



(4.2) 



where S(q' , q) is the (open) action integral along the orbit 
£ (at fixed energy e) 



S(r',r) = f p £ (r")-dr' 

J r' 



(4.3) 



projected on the q axis. The function S(q',p) is actually 
the generating function of the Poincare map: 



S(q',p)- {q,p) ->■ (q',p') ■ 



(4.4) 



If one expands the function $(t/',p) = S(q' ,p) — q'p in 
the phase of the integrand of (|4. 1|) around q' = p = 
up to quadratic terms in q' and p and evaluates the inte- 
grals in (|4.1[) by the standard stationary phase approxi- 
mation (SPA), one obtains the contribution of the orbit 
£ to Gutzwiller's trace formula (|3.1[) with the amplitude 
p.2p . Near bifurcations of the orbit £, the SPA breaks 
down and one has to include higher-order terms in the 
expansion of $(g',p). The minimum number of terms 
needed to describe a given bifurcation scenario on the 
PSS leads to the so-called "normal forms" of $(g',p), 
which depend on the type of bifurcation. Doing the in- 
tegrations in (|4.1[) using such normal forms leads to a 
finite combined contribution of the orbit £ and all the 
orbits involved in its bifurcation to the trace formula. In 
order to conform with standard notation, we rename the 
function by S(q',p) which, according to (|4.2p . is 

identical with the projected action integral S(q',q), but 
taken as a function of the variables q' and p. 
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B. Normal form for codimension-two bifurcation 

For the description of the codimension-two bifurcation 
scenario of the orbit 3 A and its satellites 3C, 3D dis- 
cussed in the previous section and illustrated on the PSS 
in Fig. I10[ the following normal form is appropriate [27] ]: 

S(I, if) = S (a) -el - al 3/2 cos(3p) - bl 2 . (4.5) 

Here one has transformed the Poincare variables {q',p) 
to quasi-polar variables (I, <p) by 



P 



21 cos ip , 



21 sin ip . 



(4.6) 



In (14. 5p . So (a) is the (closed) action integral of the central 
(3A) orbit as a function of a. The "bifurcation param- 
eter" e is a monotonously decreasing function of a such 
that e = at the touch-and-go bifurcation point (here 
a = ao = 0.407) of the central orbit, and that e < for 
a > cxq. a and b are parameters which are specific for 
the system and will be determined below. 

That (|4.5p is able to describe the correct fixed-point 
structure not only of the 3A orbit but also of its satellites 
3C and 3D, will now be shown explicitly. We first rewrite 
(|4.5p in terms of q' and p: 



> 2 



V8 X 



S(q,p) = So 
The stationary-phase conditions are 
= 0, 



4=( P 3 -3 TO 'vJ( P 2 - 



dS_ 

dq 7 



dS_ 

dp 



0. 



y 2 ) 2 . 

(4.7) 



(4.8) 



One of the solutions is qo — 0, po = and corresponds 
to the central 3 A orbit. This is so by default, due to 
the choice of the Poincare variables (q,p). Two further 
stationary points are found to be 



?i,2 = 0, 



Pi, 2 



3a 
AV2b 



± 



9a 2 
3262 



(4.9) 



The fixed points (qi,Pi) with i = 1,2 belong to the satel- 
lite orbits 3C and 3D, respectively. Two more pairs of 
fixed points for each of them are found by rotations in 
the (I,<p) plane: ip — ► ip ± 27r/3. e\ — 9a 2 /326 is the 
critical value for p 12 to have real values, i.e., for the real 
orbits 3C and 3D to exist. For our system we can choose 
b > 0, and therefore e\ > 0. For e > ex, the orbits 3C 
and 3D become imaginary and only the central orbit 3A 
is real. With decreasing e, the three pairs of stable and 
unstable satellite orbits appear at e — ex which we iden- 
tify with the tangent bifurcation point. At e = we have 
the touch-and-go bifurcation as stated above. For e < 
(a > ao), all orbits are real. 

The normal-form parameters e, a and b can be deter- 
mined by fitting the actions S and stability traces Tr M 
of the orbits, which have been obtained numerically, to 



their local behaviors predicted by the normal form (|4.5p 
The stability trace is given in terms of S by [j"c| 



TrM 



d 2 S 



dpdq 1 J 



r 



1 



d 2 S \ d 2 S d 2 S 



dpdq' J dp 2 dq' 



For the central orbit, it becomes 

TrM = Tr M A = 2 



(4.10) 



(4.11) 



so that e can be determined as e = ±\/2 — Tr Ma, choos- 
ing the correct sign on either side of the bifurcation. The 
other parameters are obtained from the action difference 
of the satellite orbits 1 and 2 (i.e., 3C and 3D). Inserting 
(qx,2,Px,2) from (|4.9p into (|4.7p . one finds 



Sd - S c 4 3/2 
^kR Q ~ = 3b Vriie - ei) ■ 



(4.12) 



By fitting the numerical data for Sd — Sc as function 
of e, we obtain e\ and b, and therefore a. Thus we can 
uniquely determine all the normal form parameters. The 
results for the triangular orbits discussed above are 



0.519252 
\fWkRo 



2.34950 

TikRo 



ei= 0.0322755. (4.13) 



The formulae for the stability traces for the orbits 3C 
and 3D are 



Tr M c ,d =2±12Jef(ex - e) - 2Ae 1 (e 1 - e) 



±12^/ex(ex-e) 3 , (4.14) 
and the action difference between orbits 3C and 3A is 
Sc — Sa 



TikRo 



1 

12b 



3e 2 - 12eie + 8e 2 - 8-vAi(ei - e) 3 



(4.15) 

We have checked that the numerical results in the neigh- 
borhood of the bifurcations are nicely reproduced by 
these equations. 



C. Uniform approximations 



Inserting the normal form (|4. 5[) into the integral (|4.ip , 
one obtains a "local" uniform approximation 8] which 
is finite and valid near the bifurcation, i.e., for not too 
large absolute values of e. Due to the C% v symmetry of 
our system, the touch-and-go bifurcation is non-generic 
and isochronous. The factor n r in the denominator of 
(|4.ip here must be chosen [28| as n r — f = 3. We have 
another degeneracy factor 2 in the numerator due to the 
time-reversal symmetry of all orbits. Using the variables 
(/, ip) and changing the energy e to the wave number fc, 
we finally obtain the following expression for the local 
uniform approximation, in which the integration over ip 
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can be done analytically: 

s ^ = ^n Re l dI J d ^ 



d 2 s 



1/2 



dldip 

( So -lip -el- al 3/2 cos(3<^) - W 2 } 
u± Re e lkL t I 2 / d/Jo ( |/ 3/2 )e 



x exp 



-b/ 2 ) 



(4.16) 



The integration over / can be performed numerically us- 
ing an expansion formula given in [27| . 

As stated above, the result (|4.16p is only useful in 
the neighborhood of the bifurcation. Far away from it, 
where all involved orbits become isolated, it can be eval- 
uated asymptotically (corresponding to the SPA), but 
the amplitudes of the orbits then do not agree with the 
Gutzwiller values (|3 . 2|) . In order to achieve this, one 
must develop "global" uniform approximations [H, Hc| |. 
To that purpose, one needs to include higher-order ex- 
pansion terms in the normal form. For the codimension- 
two bifurcation of our type one obtains, after suitable 
coordinate transformations [2al , 



%(fc) 



x exp 



with 



{So 



Re 



dl 



d<p* (I, ip) 



lip - el 



a / 3 / 2 COs(3(^)-W 2 }-y^ 



<p) = L 6 -aI- f3I 3/2 cos(3cp), 



(4.17) 



(4.18) 



where a and (3 are expressed by a certain combination 
of higher-order expansion coefficients. In practice, these 
parameters are determined such that (|4.17p yields the 
sum of contributions with amplitudes (|3.2p of all involved 
isolated orbits far away from the bifurcation point. They 
are determined by equating (cf. [30() 



L 



c 



La - ah 



^|2-TrM c | 
Ld 



La - alo 



y/\2-TrM D \ vTdet $"(I D )\ 



with 



and 





d 2 S 


d 2 S 




d 2 S 


d 2 S 


det = 


dip' 2 
d 2 S 


dipdl 
d 2 S 






dq' dp 
d 2 S 




OWip 




ip=0 


dpdq' 


dp 2 



q>=0 



'C,D 



<?1,2 



(4.19) 
(4.20) 

(4.21) 

(4.22) 



We now have all ingredients ready for calculating the 
semiclassical level density of our system. Since it is well- 
known l| that the sum over all periodic orbits does not 
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FIG. 11: Same as Fig. [5]for the hexagonal orbits. 

converge in systems with mixed dynamics, we have to 
introduce a certain truncation. This is achieved @,[3l[ by 
focusing on the gross-shell structure of the level density, 
coarse-graining it by a Gaussian convolution over k with 
width 7: 

9l {k) = -J=-Jdk' e^ k '- k ^ 2 g(k') . (4.23) 
The semiclassical trace formula (13.31) then becomes 



9j(k) ^ g(k) + XX 



h L i/ 2 ) 2 A^(k)cos(kL 6 - - Vi ) 



(4.24) 

Due to the Gaussian factor, periodic orbits with lengths 
L^ > 27r/7 are now exponentially suppressed. Choos- 
ing 7 = 0.6, we need only consider periodic orbits with 
< 10i?o- From Fig. [H we see that no isochronous 
bifurcations of the diametric orbit occur for a > 0, so 
that its contribution can be evaluated by the standard 
Gutzwiller formula with amplitudes (|3.2[) . A period- 
doubling bifurcation of the stable diameter occurs at 
a ~ 0.15, but does not contribute much to the coarse- 
grained level density. The same is true for the tetragonal 
and pentagonal orbits. The contributions of the trian- 
gular orbit 3 A and its satellites 3C, 3D are included in 
the global uniform approximation (|4.17[) . including the 
Gaussian damping factor exp[— (^/La/2) 2 ]. 

The hexagonal orbit encounters also a codimcnsion- 
two bifurcation of the same type as the triangular orbit, 
as seen in Fig. 1111 Determining the normal-form param- 
eters in the same manner as above, we find the values 



14.0046 



414.669 

flkRn 



y/nkRo nkR 

which are much larger than those for the triangular orbit 
given in (|4.13p . Since they contribute inversely to the 
local uniform approximation (|4.16[) . this bifurcation has 
a much less dramatic influence on the level density than 
that of the triangular orbit. In order to demonstrate this 
more explicitly, we plot in Figs. [T^] and fTS] the modulus 
of the integral 



AOC 

F(e;a,b)= / dIJ 
Jo 



3/2x i(-d-bl 2 ) 



(4.26) 




b 



—e — e 



FIG. 12: Profile of the function F(e; a, b) in for the 

normal form parameters (|4.13[) of the triangular orbit. The 
left panel is for fixed b and varying a = a^/hkRo, and the right 
panel for fixed a and varying b — bhkRo, both for kRo = 20. 
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FIG. 13: Same as Fig. 1121 but with the normal form param- 
eters (|4.25[1 of the hexagonal orbit. 



appearing in (|4. 16|) . Fig. [12] is for the parameters (|4.13[) 
of the triangular and Fig. [T3] for the parameters (|4.25p 
of the hexagonal orbits. We see that for larger a and b, 
this integral has smaller values and is a more monotonous 
function of e. For small a and b, however, it takes large 
values and exhibits a considerable peak near the bifur- 
cation point e = 0. (Note that the maximum actually 
occurs slightly after the bifurcation, i.e., for — e > 0, as 
discussed at the end of Sect. Hill ) This can be under- 
stood as follows. If the normal-form parameters a and b 
are small, the action S does not depend much on I and ip, 
and the quasiperiodic orbits around the central orbit give 
a more coherent contribution to the integral, yielding a 
better restoration of local dynamical symmetry around 
the central periodic orbit. Actually, the limit a — » is 
sufficient to locally restore integrability: the normal form 
(|4.5p then only depends on I and for 6^0 describes the 
bifurcation of a torus from an isolated orbit (cf. 32]). 

This local dynamical symmetry is equivalent to a res- 
onance condition for quasi-tori winding around the bi- 
furcating periodic orbit over a rather wide region. The 
quantization of these quasi-tori yields a number of quasi- 
degenerate levels, seen as the bunches in the spectrum 
of Fig. [3] This increases the probability for small level 
spacings and renders the NNS distribution more Poisson 
like. The considerable changes in the NNS distributions 
shown in Figs. [5] and [5] around a = 0.5 therefore account 
for the emergence of these quasi-resonant tori. 

In Fig. [14] we compare the oscillating part of the 
semiclassical level density with that of the quantum- 
mechanical result. In the semiclassical calculation we 
have included the primitive diametric, triangular, tetr ag- 
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FIG. 14: Oscillating part of the coarse-grained level density. 
Solid and dashed lines represent the semiclassical and the 
quantum-mechanical results, respectively. 



onal, pentagonal and hexagonal orbits as explained 
above. We see that, besides the rapid oscillations due to 
the average length of the shortest orbits included (i.e., 
mainly of the short diameter orbit), there is a beat- 
ing pattern that comes from the interferences between 
the different orbits. We note that both amplitudes and 
phases are nicely reproduced in the semiclassical result 
for all deformation parameters including the bifurcation 
points a ~ 0.4. For a = 0.5 and 0.6, the interference 
effects are reduced and the oscillating pattern is quite 
regular, indicating the dominance of the bifurcating pe- 
riodic orbits and the accompanying quasi-tori which have 
approximately the same lengths. This is also seen in the 
shell-correction energy for a = 0.5 in Fig. [U where the 
spacing AN of the regular oscillation corresponds to the 
length of the dominating triangular orbit. 



V. DISCUSSION OF DISCRETE SYMMETRIES 

We conclude by some remarks on the role of discrete 
symmetries. In the above model, the triangular periodic 
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orbit undergoes a bifurcation in the course of which a 
pair of three-fold degenerate new orbits are created. (We 
only mention the relative degeneracies here; all orbits 
have two extra degeneracy factors of 2 due to reflections 
at the x axis and to time-reversal.) Stated more gener- 
ally: the bifurcation brings about new periodic orbits of 
reduced symmetry, which have several degenerate repli- 
cas, connected with the symmetry operations 1Z and V in 
(|2.4p . that give coherent contributions to the level den- 
sity. This is one of the reason why we obtain a strong 
shell effect due to this bifurcation. If the C3 symme- 
try is slightly broken, the equilateral triangular orbit will 
split into 3 non-equilateral triangular orbits with differ- 
ent lengths and different stabilities. They will give de- 
structive contributions to the level density in most of the 
parameter region. We would therefore expect that in bil- 
liards with C nv symmetry, the regular polygonal orbits 
with n reflections will play the most prominent role, such 
as the triangular orbit in the present system. 

This behavior can also be predicted from a perturba- 
tive trace formula, developed by Creagh [33j], which de- 
scribes semiclassically the breaking of continuous sym- 
metries. In systems with continuous symmetries, the 
leading periodic orbits occur in degenerate families cor- 
responding to rational tori. When a perturbation breaks 
the continuous symmetries, the tori are broken into iso- 
lated orbits. For weak perturbations, only those orbit 
families (tori) £ are broken for which the action change 
A5j in lowest order of the classical perturbation theory 
is nonzero. Furthermore, if the perturbed system still 
has a discrete point symmetry, the orbit families which 
have the same point symmetry are in resonance with the 
perturbation and typically suffer the largest first-order 
action change. In our billiard system (|2.2[) . we can treat 
the deviation from the circular billiard, for sufficiently 
small values of a, as a perturbation that breaks the con- 
tinuous U(l) symmetry. For small a, the shape of the 
boundary is given by 

R(9) «-Ro(l-£cos30), (5.1) 

and e = v3c*/9 is the appropriate perturbation param- 
eter. As we show in Appendix [B] the length of the pe- 
riodic orbit family (v, w) is changed in first order of e 
only for v = 3w, i.e., only triangular orbits are affected 
in first-order perturbation theory. The situation is simi- 
lar in a spherical cavity perturbed by deformations with 
2 J -pole deformations. For this system it was shown ex- 
plicitly [iH that the orbit families with a j-fold sym- 
metry obtain the largest first-order action changes. As 
a consequence of this symmetry-breaking argument, we 
understand that the orbits with the highest point sym- 
metry (namely that of the perturbed system itself) are 
those whose stability traces TrMj deviate fastest from 
their values Tr = +2 at a — 0. Their stable branches 
therefore are also the fastest to undergo a bifurcation. 

We finally note that, quite generally, orbits with point 
symmetries keep undergoing bifurcations whereby the 
new-born orbits at each successive bifurcation lose one 
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FIG. 15: Trace of stability matrix for the square orbits (up- 
per panel) and the second repetition of the diametric orbits 
(lower panel) in the circular-to-square billiard (|5.2p with Cav 
symmetry. 

symmetry, until the last stable branch has lost all sym- 
metries. (We have no formal proof for this statement, but 
have observed this scenario in many Hamiltonian sys- 
tems.) Hereby the stable branches survive up to large 
deformations, giving large contributions to the level den- 
sity, while the unstable branches rapidly turn chaotic and 
their contribution decreases exponentially. In conclusion, 
orbits with high point symmetries tend to give large con- 
tributions to the level density, even in well-deformed sys- 
tems, as long as they possess discrete symmetries. 

To illustrate some of the above statements, let us con- 
sider a billiard system with a different discrete symmetry. 
The square billiard is integrable, and we can smoothly 
connect square and circular billiards by parameterizing 
the boundary shape R{9) as 

„ 2 ai? 4 (l + cos46>) - Rk ~ . , 



For small a, R(9) is expressed by R{6) ~ i? (l + ecos4#) 
with e = a/16. This system obviously has C^ v symmetry. 
a = corresponds to the circle with radius Rq , and a = 1 
corresponds to the square with side length -JT/^Rq. The 
system is again integrable at these two limits, while it is 
non-integrable in between. 

The square orbit here also has C iv symmetry, and it 
suffers a non-generic island-chain bifurcation at a = 0.5, 
as seen in the upper panel of Fig. [T5] (see Appendix [A] 
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FIG. 16: Same as Fig. [7] for the circular-to-square billiard 
(|5.2|) with d v symmetry. 



for the analytical form of its stability matrix) . The pair 
of new-born stable and unstable orbits have lost the C^ v 
symmetry but still have a reflection symmetry at one 
axis, they are therefore doubly degenerate. The stable 
branch suffers a pitchfork bifurcation at a ~ 0.8, where 
it loses the last reflection symmetry. The diameter orbit 
(seen in the lower panel of Fig. [15]) has Ci v symmetry; 
its second repetition suffers an island-chain bifurcation at 
a = 0.5, where a new pair of stable and unstable branches 
are born. At this point, the stable branch loses its time- 
reversal symmetry and keeps the Civ symmetry. The or- 
bit born at a pitchfork bifurcation near a ~ 0.7 loses the 
reflection symmetry but still keeps C2 symmetry. It suf- 
fers a further pitchfork bifurcation near a ~ 0.85 and the 
orbit born there has lost all symmetries. The octagonal 
orbits also show interesting bifurcation sequences in this 
system, but their contributions to the level density are 
small compared to those of the above ones. Other short 
orbits without symmetries are already strongly chaotic 
at much smaller values of a. 

Figure [16] shows the modulus of the Fourier trans- 
form (|3.6p of the quantum-mechanical level density of 
the circular-to-square billiard (|5.2[) . The above bifurca- 
tions are seen to play a significant role, although their 
effects are not so drastic compared as for the circle-to- 
triangle billiard, since the square and diameter orbits give 
destructive contributions. 



VI. SUMMARY AND CONCLUSIONS 

We have used periodic orbit theory to analyze a strong 
enhancement of shell effects in a non-integrable billiard 
system, which is closely connected to the bifurcation of 
a periodic orbit with high point symmetry. The semi- 
classical trace formula, including the bifurcating orbits 
in a global uniform approximation and the other isolated 
orbits with their standard semiclassical Gutzwiller am- 



plitudes, describes very well the coarse-grained quantum 
level density of the system. We have shown that the 
bifurcations play a significant role especially when their 
normal-form parameters are small. Under such a con- 
dition the quasi-periodic orbits around the bifurcating 
orbit give coherent contributions to the trace integral, 
which can be considered as an approximate local sym- 
metry restoration around the periodic orbit. It is accom- 
panied by a relatively large regular region in the phase 
space which persists around a ~ 0.5 (cf. Fig. [TU]) and 
is reflected by the level statistics (Figs. [5] 0, and by 
large shell effects in the spectrum (Fig. [3]) and the shell- 
correction energy (Fig. [4]). 

The role of the discrete symmetry is also important for 
creating a large shell effect. The bifurcation of a highly 
symmetric orbit causes degenerate symmetry-reduced or- 
bits which give additional contributions to the level den- 
sity. Through successive bifurcations, at each of which 
one symmetry is lost by the new-born orbits, their sta- 
ble branches may survive up to rather large deformations 
until all symmetries are lost and all remaining branches 
become unstable with exponentially decreasing contribu- 
tions to the level density. This mechanism is operative 
in systems with arbitrary point symmetries (such as C nv 
discussed here) , in which periodic orbits of the same sym- 
metries exist. 
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APPENDIX A: ANALYTIC EXPRESSIONS FOR 
TRACES OF STABILITY MATRICES 

If the geometry of a periodic orbit in a billiard is an- 
alytically given, analytic expressions for its stability ma- 
trix are often available as well. In billiard problems, the 
stability matrix is constructed by alternating products of 
translation and reflection matrices (see, e.g., Appendix B 
of [24j). For a periodic orbit with v vertices, it becomes 

M = M Te{ {v)M tTa (v, v -!)■■■ M rof (l)Af tra (l, 0) , (Al) 

where M re f(i) stands for the matrix for reflection at the 
i-th vertex, and Mt Ta (i + X,i) that for translation from 
the i-th to the i + 1-th vertex, respectively. They are 
given analytically by 

M tia (i + l,i)=(l Li+ Y ,P ^j, (A2) 

^0 = (v*-<* -1)' (A3) 
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where Lj+i,, is the distance between the vertices i and 
p the momentum, pi the curvature radius of the 
boundary at i-th vertex, and <fi the reflection angle mea- 
sured from the normal to the boundary at the i-th vertex. 
For a boundary R(9), the curvature radius is given by 



P = 



R 2 + R" 



3/2 



R 2 + 2R' 2 - RR" 



(A4) 



where primes indicate derivatives with respect to 9. 

We first apply this to the equilateral triangular orbit 
3 A in the billiard system (|2.2[) . Its vertices are located at 
8 = 0, ±27r/3, and the curvature radius at these points is 



(3 - 7 2 b 
2(4 7 2 - 3) 



Rn 



(A5) 



where 7 is the distance between the origin and the reflec- 
tion points in units of Rq] it is a function of the defor- 
mation parameter a given implicitly by the equation 



7 2 + 



2V3a 



9 



-7 



= 1 . 



(A6) 



Since the reflection and translation matrices for all ver- 
tices and sides are equal, the stability matrix is simply 
given by M = (M ref M tr ) 3 = Mf. Using the symplectic 
nature of the matrices M, we can express the trace of the 
full M in terms of the trace of M\: 



TrM = Tr(A/i) 3 = (Tr Mi) 3 - 3Tr Mi . 



(A7) 



Inserting L i+ is = ^/3jRq, fi = tt/6 and (|A5|) into (|A2|) 
and (EH, we'get TrMi = (34 7 2 - 30)/(3 - j 2 ). At 
the bifurcation point, TrMi = — 1 and thus 7 = 3/vTT. 
Inserting this into (|A6I) . we get the deformation 



a 3A = 11/27= 0.407407. 



(A8) 



Bifurcations of the square orbit 4A and the linear or- 
bit 2A in the circle-to-square billiard system (|5.2p are 
analyzed in the same way. For both orbits, we have 
TrMi = 2 -4a. Thus, orbit 4A and the 2nd repetition of 
orbit 2 A have the same value of Tr M at all deformations. 
The bifurcation deformation becomes, from TrMi = 0, 



C*4A = «2A 2 = I/ 2 ■ 



(A9) 



APPENDIX B: CLASSICAL PERTURBATION OF 
THE CIRCULAR BILLIARD 



by a deformation 



R{9) =Ro(l-e cos n6 



(Bl) 



We consider the primitive periodic orbit family (v,w), 
with mutually prime integers v and w representing the 
number of vertices and the winding number around the 
origin, respectively. Due to the deformation, the posi- 
tions of the vertices 9i (i = 1, ■ ■ ■ ,v) are shifted by 56i, 
which are of order e. The distance between successive 
vertices is given by 



Lf+i i = R\ + Ri+i — 2RiRi+i cos(8 i+1 — 9i) , 



Ri — R(9i), 9i — 9q + i 



2W7T 



(B2) 



where one of the two branches has 9q = and the other 
branch has 9q = ir/l nv , l nv being L.C.M. of n and v. 
Inserting (|B1[) and keeping only terms up to first order 
in e, we obtain 



W7T 



AL !+ i, w R ^ cos — (59 i+ i - 59i) 

2niunr\ 



v 

WTT 



e sin ■ 



v 



cos n9, 



70 



— 



cos n9 { 



2n(i + l)wn 



(B3) 



Therefore, the changes of the periodic orbit lengths are 



AL VW — S ' A_Li+i : 



-2eR sin — 2^ cos ( n9 H ) . (B4) 

2—1 



The sum on the right-hand side has a non- vanishing value 
only for n/v = t,t being an integer; then it becomes 



WTT WTT 

AL VW = —2vsRq sin cos{n9 Q ) = =f2veRo sin . 

v v 

(B5) 



In this appendix we calculate the change of the periodic 
orbit lengths in the circular billiard, when it is perturbed 



All other periodic orbit families suffer a change of length 
only at higher orders of s. 
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